# include "scene.h"
// # include "discr.h"
# include <stdio.h>
# include <stdlib.h>

#ifdef _OMP_PRESS_BC_
# include <omp.h>
#endif

extern double inlet_velocity(double r, double t, short index);
// extern short interpol_1dim(double const *x, double const *u, unsigned short n, double xp, double *up);
/*
# ifdef DENS_TDEP
extern double dens(double T);
# endif
# ifdef KVISC_TDEP
extern double kvisc(double T);
# endif*/
/*
# include <time.h>
# include <sys/timeb.h> */

extern void press_bc(struct grid *g, short tflag, unsigned short iRK)
  {/* tflag = 0: current time step i.e. t; tflag = 1: next time-step i.e. t+dt. */
   short i, j, k, nr, nz, nth, /*tid,*/ nTRD, indx;
   double ***u, ***v, ***w, ***p,/* ***T,*/ *r, *th, *z, /*L, */density0, t, /*kvisc0,*/ /*_z[PRESS_INTRP_N], **_p,*/ q1, coeff;
   struct cmtrx *cm;
// struct timeb t1, t2;

// ftime(&t1);

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 2) )
       {
        fprintf(g->lptr, "\nERROR: In 'press_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

   density0 = g->density0; /* kvisc0 = g->kvisc0; */ nth = g->nth;  th = g->th;
   t = g->t + tflag*g->dt;

   nTRD = 1;
#  ifdef _OMP_PRESS_BC_
   nTRD = g->nTRD;
#  endif

/* These are used in the interpolation of pressure at outlet. */
//    _p = (double **)malloc(nTRD*sizeof(double *));
//    for(i = 0; i < nTRD; i++) _p[i] = (double *)malloc(PRESS_INTRP_N*sizeof(double));

/* ATTENTION: T is utilized to calculate physical properties of the fluid (kvisc, density etc.) wherever Neumann-BC
   for the pressure is applied during course of the calculation of Navier-Stokes momentum terms. */

/* ----------------- Applying boundary condition for SECTION - 0 ------------------------- */
   u = g->u[0]; v = g->v[0]; w = g->w[0]; p = g->p[0]; /*T = g->T[0];*/ r = g->r[0]; z = g->z[0]; nr = g->nr[0]; nz = g->nz[0];

#  ifdef _OMP_PRESS_BC_
   #pragma omp parallel private(i,j,k,/*L,tid*/,q1,coeff,cm)
#  endif
     {
   /* At the top boundary open to the section-1. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++)
           for(k = 0; k < NzG; k++) p[i][j][NzG+nz+k] = g->p[1][i][j][NzG+k];

   /* At the INLET: Numan BC */
//       k = NzG-1; /*** For the BOTTOM INTERFACE. It is used in evaluation of "L" while discretizing manually. ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr_nozzle; i++) for(j = NthG; j < nth + NthG; j++)
           {
         /* Calculating [dw/dt] for iRK-th Runge-Kutta step. */
            switch(iRK)
               {
                case 0:
                   q1 = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1);
                   break;

                case 1:
                   q1 = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                         + c2*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2);
                   break;

                case 2:
                   q1 = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                         + (a31+a32)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2)\
                         + a32*c2*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 3);
                   break;
               }

//               L =  0 - q1 - ZMOMT_TERM_3 + ZMOMT_VISC + ZMOMT_FORCE; /* Doesn't contain any "p" term. */

            q1 = density0*(g->C[0][i-NrG][j-NthG][0] - q1);

            cm = &g->dcmx_zcs1[0][0];
            for(indx = 0; indx <= cm->iend[NzG-cm->respadi]-cm->istart[NzG-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NzG-cm->respadi]+indx == NzG-1) coeff = cm->mtrx[NzG-cm->respadi][indx];
               else q1 -= cm->mtrx[NzG-cm->respadi][indx]*p[i][j][cm->srcpadi+cm->istart[NzG-cm->respadi]+indx];

            p[i][j][NzG-1] = q1/coeff;

//             p[i][j][NzG-1] = p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
           }

   /* At the SOLID-BASE after the nozzle inlet: Numan BC */
//       k = NzG-1; /*** For the BOTTOM INTERFACE. It is used in evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_zcs1[0][0];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NzG-cm->respadi]-cm->istart[NzG-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NzG-cm->respadi]+indx == NzG-1) coeff = cm->mtrx[NzG-cm->respadi][indx];
               else q1 += cm->mtrx[NzG-cm->respadi][indx]*p[i][j][cm->srcpadi+cm->istart[NzG-cm->respadi]+indx];

//              L = ZMOMT_VISC + ZMOMT_FORCE; /* Doesn't contain any "p" term. */

             p[i][j][NzG-1] = (density0*g->C[0][i-NrG][j-NthG][0]-q1)/coeff;

//              p[i][j][NzG-1] = p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
            }

   /* At the cylindrical SIDE-WALL: Numan BC */
//       i = NrG+nr-1; /*** For the cylindrical-wall. It is used in evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_rcs1[0];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NrG+nr-cm->respadi]-cm->istart[NrG+nr-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx == NrG+nr) coeff = cm->mtrx[NrG+nr-cm->respadi][indx];
               else q1 += cm->mtrx[NrG+nr-cm->respadi][indx]*p[cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx][j][k];

             p[NrG+nr][j][k] = (density0*g->A[0][nr][j-NthG][k-NzG]-q1)/coeff;

//              L = RMOMT_VISC; /* Doesn't contain any "p" term. */

//              p[NrG+nr][j][k] = p[NrG+nr-1][j][k] + density0*(r[NrG+nr]-r[NrG+nr-1])*L;
            }

   /* For the NEGATIVE radius: applying symmetric B.C. across the center. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(dynamic,3)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
          p[NrG-1-i][j][k] = p[NrG+i][j][k];

  /* For ghost-cells in theta-direction: Applying Periodic BC. */
#    ifdef _OMP_PRESS_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         p[i][NthG-1-j][k] = p[i][NthG+nth-1-j][k];
         p[i][NthG+nth+j][k] = p[i][NthG+j][k];
       }


/* ----------------- Applying boundary condition for SECTION - 1 ------------------------- */
#     ifdef _OMP_PRESS_BC_
      #pragma omp single
#     endif
        {
         u = g->u[1]; v = g->v[1]; w = g->w[1]; p = g->p[1]; /*T = g->T[1]; */r = g->r[1]; z = g->z[1]; nr = g->nr[1]; nz = g->nz[1];
        }

   /* At the lower interface open to the section-0. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < NrG+g->nr[0]; i++) for(j = NthG; j < NthG+nth; j++)
         for(k = 1; k <= NzG; k++) p[i][j][NzG-k] = g->p[0][i][j][NzG+g->nz[0]-k];


   /* At the SOLID-BASE between the interface-open-to-the-lower-section-0 and the outlet: Numan BC. */
//       k = NzG-1; /*** For the BOTTOM INTERFACE. It is used in evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr[0]; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_zcs1[1][1];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NzG-cm->respadi]-cm->istart[NzG-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NzG-cm->respadi]+indx == NzG-1) coeff = cm->mtrx[NzG-cm->respadi][indx];
               else q1 += cm->mtrx[NzG-cm->respadi][indx]*p[i][j][cm->srcpadi+cm->istart[NzG-cm->respadi]+indx];

//              L = ZMOMT_VISC + ZMOMT_FORCE; /* Doesn't contain any "p" term. */

             p[i][j][NzG-1] = (density0*g->C[1][i-NrG][j-NthG][0]-q1)/coeff;

//              p[i][j][NzG-1] = p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
            }

   /* At the outlet open to the extended section-2. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+nr-g->nr_gap; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
          for(k = 1; k <= NzG; k++) p[i][j][NzG-k] = g->p[2][i-nr+g->nr_gap][j][NzG+g->nz[2]-k];

   /* At the cylindrical SIDE-WALL: Numan BC. */
//       i = NrG+nr-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_rcs1[1];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NrG+nr-cm->respadi]-cm->istart[NrG+nr-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx == NrG+nr) coeff = cm->mtrx[NrG+nr-cm->respadi][indx];
               else q1 += cm->mtrx[NrG+nr-cm->respadi][indx]*p[cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx][j][k];

//              L = RMOMT_VISC;  /* Doesn't contain any "p" term. */

             p[NrG+nr][j][k] = (density0*g->A[1][nr][j-NthG][k-NzG]-q1)/coeff;

//              p[NrG+nr][j][k] = p[NrG+nr-1][j][k] + density0*(r[NrG+nr]-r[NrG+nr-1])*L;
            }

   /* At the TOP-SURFACE: Numan BC. */
//       k = NzG+nz-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(dynamic,1)
#     endif
      for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_zcs1[1][0];
             if( (i >= NrG+g->nr[0]) && (i < NrG+g->nr[1]-g->nr[2])) cm = &g->dcmx_zcs1[1][1];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NzG+nz-cm->respadi]-cm->istart[NzG+nz-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NzG+nz-cm->respadi]+indx == NzG+nz) coeff = cm->mtrx[NzG+nz-cm->respadi][indx];
               else q1 += cm->mtrx[NzG+nz-cm->respadi][indx]*p[i][j][cm->srcpadi+cm->istart[NzG+nz-cm->respadi]+indx];

//              L = ZMOMT_VISC + ZMOMT_FORCE; /* Doesn't contain any "p" term. */

             p[i][j][NzG+nz] = (density0*g->C[1][i-NrG][j-NthG][nz]-q1)/coeff;

//              p[i][j][NzG+nz] = p[i][j][NzG+nz-1] + density0*(z[NzG+nz]-z[NzG+nz-1])*L;
            }

   /* For the NEGATIVE radius: applying symmetric-BC across the center. */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(dynamic,5)
#     endif
      for(k = 0; k < nz+2*NzG; k++) for(j = NthG; j < NthG+nth; j++) for(i = 0; i < NrG; i++)
          p[NrG-1-i][j][k] = p[NrG+i][j][k];

  /* For ghost-cells in theta-direction: Applying Periodic-BC. */
#    ifdef _OMP_PRESS_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         p[i][NthG-1-j][k] = p[i][NthG+nth-1-j][k];
         p[i][NthG+nth+j][k] = p[i][NthG+j][k];
        }


  /* ----------------- Applying boundary condition for SECTION - 2 ------------------------- */
#     ifdef _OMP_PRESS_BC_
      #pragma omp single
#     endif
        {
         u = g->u[2]; v = g->v[2]; w = g->w[2]; p = g->p[2]; /*T = g->T[2]; */r = g->r[2]; z = g->z[2]; nr = g->nr[2]; nz = g->nz[2];
//          for(_z[0] = z[NzG], k = 1; k < PRESS_INTRP_N; k++) _z[k] = 0.5*(z[NzG-1+k]+z[NzG+k]);
        }

#     ifdef _OMP_PRESS_BC_
//       tid = omp_get_thread_num();
      #pragma omp for schedule(guided) nowait
// #     else
//       tid = 0;
#     endif
      for(i = NrG; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
         {
       /* The top-interface open to the section-1. */
          p[i][j][NzG+nz] = g->p[1][i+g->nr[1]-g->nr_gap][j][NzG];

       /* The bottom outlet. Interpolating the pressure at the ghost-cell from the inner-cells. */
//           _p[tid][0] = P_REF + density0*(-1*GRAVITY_DOWNWARD)*z[NzG];
          cm = &g->icmx_zcs[2][0];
          q1 = 0.0;
          for(indx = 0; indx <= cm->iend[NzG-cm->respadi]-cm->istart[NzG-cm->respadi]; indx++)
             if(cm->srcpadi+cm->istart[NzG-cm->respadi]+indx == NzG-1) coeff = cm->mtrx[NzG-cm->respadi][indx];
             else q1 += cm->mtrx[NzG-cm->respadi][indx]*p[i][j][cm->srcpadi+cm->istart[NzG-cm->respadi]+indx];

          p[i][j][NzG-1] = (P_REF + density0*(-1*GRAVITY_DOWNWARD)*z[NzG] - q1)/coeff;
/*        for(k = 1; k < PRESS_INTRP_N; k++) _p[tid][k] = p[i][j][NzG-1+k];
          interpol_1dim(_z, _p[tid], PRESS_INTRP_N, 0.5*(z[NzG-1]+z[NzG]), &p[i][j][NzG-1]);*/
         }

   /* At the outer solid-wall: Neumann BC. */
//       i = NrG+nr-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_rcs1[2];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NrG+nr-cm->respadi]-cm->istart[NrG+nr-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx == NrG+nr) coeff = cm->mtrx[NrG+nr-cm->respadi][indx];
               else q1 += cm->mtrx[NrG+nr-cm->respadi][indx]*p[cm->srcpadi+cm->istart[NrG+nr-cm->respadi]+indx][j][k];

//              L = RMOMT_VISC;  /* Doesn't contain any "p" term. */

             p[NrG+nr][j][k] = (density0*g->A[2][nr][j-NthG][k-NzG]-q1)/coeff;

//              p[NrG+nr][j][k] = p[NrG+nr-1][j][k] + density0*(r[NrG+nr]-r[NrG+nr-1])*L;
            }

   /* At the inner solid-wall: Neumann BC. */
//       i = NrG-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             cm = &g->dcmx_rcs1[2];
             q1 = 0.0;
             for(indx = 0; indx <= cm->iend[NrG-cm->respadi]-cm->istart[NrG-cm->respadi]; indx++)
               if(cm->srcpadi+cm->istart[NrG-cm->respadi]+indx == NrG-1) coeff = cm->mtrx[NrG-cm->respadi][indx];
               else q1 += cm->mtrx[NrG-cm->respadi][indx]*p[cm->srcpadi+cm->istart[NrG-cm->respadi]+indx][j][k];

//              L = RMOMT_VISC;  /* Doesn't contain any "p" term. */

             p[NrG-1][j][k] = (density0*g->A[2][0][j-NthG][k-NzG]-q1)/coeff;

//              p[NrG-1][j][k] = p[NrG][j][k] - density0*(r[NrG]-r[NrG-1])*L;
            }

  /* For ghost-cells in theta-direction: Applying Periodic BC. */
#    ifdef _OMP_PRESS_BC_
     #pragma omp for schedule(dynamic,3)
#    endif
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
         p[i][NthG-1-j][k] = p[i][NthG+nth-1-j][k];
         p[i][NthG+nth+j][k] = p[i][NthG+j][k];
        }
     }

//    for(i = 0; i < nTRD; i++) free(_p[i]); free(_p);

   return;
  }
